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■ Abstract 

CO , The motion of binary star systems is re-examined in the presence of 

^ ' perturbations from the theory of general relativity. To handle the singu- 

larity of the Kepler problem, the equation of motion is regularized and 
' , linearized with quaternions. In this way first order perturbation results 

' are derived using the quaternion based approach. 

^ ■ 

O ■ 1 Introduction 

(N ; 

" , I . In this paper gravitational effects as perturbations of the Kepler problem are 

^ ' examined with post-Newtonian approximation. Gravitational effects become 

. strong when the components of the binary are close to each other, and the orbital 

' separation is sman0 The Kepler problem is singular when the separation is zero, 

. therefore to study gravitational effects the desingularization - or regularization 

- of the equation of motion would be a substantial step. 

It is well known that Kustaanheimo (1964) solved the regularization of the 
three-dimensional Kepler problem with spinors[T], which was reformulated by 
Stiefel[2]. In their method - the KS method for short - the regularization was 
carried out in four dimensions, and it was proved that the three-dimensional 
Kepler problem can only be regularized using four-dimensional linear spaces. 

We follow another approach developed by J. Vrbik. In his work the men- 
tioned four-dimensional space is the linear space of quaternions and the regu- 
larization is calculated with quaternion algebra. He applied his method with 



^The post-Newtonian approximation is not valid when the orbital separation is smaller 
than the innermost stable circular orbit or when the bodies start to merge. Therefore in the 
paper it is supposed that these limits are not reached. 
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success to the Lunar probleni[5], and several perturbative forces were studied in 
details H [316]. 

In the present work we use his method to examine gravitational effects an- 
alytically with quaternions. The leading order correction of general relativity 
to classical mechanics is calculated first. The formula for the precession of 
the pericentre is derived based on the Vrbik's quaternion formulae. Then the 
gravitational radiation reaction is analyzed, where the famous Peters-Mathews 
formula is proved[7]. In this calculation we manage to remove the residual 
coordinate gauge freedom of the gravitational reaction from the quaternionic 
equation of motion. In addition using a one-dimensional model we demon- 
strate that the regularization can lead to different results depending on that 
the Sundman transformation is employed with the perturbed or unperturbed 
orbital separation. 

The regularization is defined with four-dimensional spaces, thus an addi- 
tional geometrical constraint - a gauge - have to be applied to describe the 
three-dimensional spatial Kepler problem. In the KS method the so-called bi- 
linear relation is defined, which is an excellent gauge for numerical calculations. 
Vrbik proposes another constraint to provide an analytic perturbative method, 
since - according to Vrbik - the bilinear relation is too restrictive to build an 
analytic perturbative method. This constraint is the major difference between 
the KS method and Vrbik's work. 

The Laplace vector is a constant of motion of the Kepler problem, which is a 
consequence of the hidden symmetry of the problem[8]. This symmetry becomes 
manifest in four dimensions, which shows that the Kepler problem has another 
interesting connection with the four-dimensional space. This connection has far 
reaching consequences P ITO] . 

Quaternions were first applied to regularize the Kepler problem by Chel- 
nokov who successfully regularized the Kepler problem to describe rotating co- 
ordinate systems [11]. Moreover he was able to apply his results to describe the 
optimal control problem of a spacecraft [12] . 

Later it was shown by Vivarelli (1983) in a general mathematical sense 
that the KS method can be transformed identically into quaternion algebra [T5]. 
Quaternion algebra proved to be very useful to derive the central ideas of the 
KS method. Remarkably the bilinear relation is described as a fibration of the 
quaternion space. 

More recently Waldvogel showed that the spatial Kepler motion can be ele- 
gantly formulated with quaternions using a novel star conjugation operator|14j. 
The star conjugation is especially useful to handle the bilinear relation. The in- 
teresting connection with the BirkhofF transformation is also shown [15] . Quater- 
nions turned to be useful in case of three and N-body applications [TB] . 

It has to be emphasized that the mentioned quaternion approaches exclu- 
sively apply the "bilinear relation" as a gauge to reduce the dimensions from 
four to three, while Vrbik apply his special gauge. 

This paper is organized as follows: a short outline of Vrbik's approach is 
provided in Section [2] and [S] where we describe the transformation of the Kepler 
problem into quaternion differential equation. Then the solution is given in 
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terms of ordinary differential equations of orbital elements. The advantages of 
Vrbik's calculus compared with the KS method are highlighted. 

In Section r2.5l a onc-dimcnsional example is given where we demonstrate that 
the result of the regularization depends on whether the Sundman transformation 
is applied with perturbed or unperturbed orbital separation. 

In Section!?] Vrbik's method is applied to two perturbations. First of all, the 
leading order correction of general relativity to classical mechanics is examined. 
The formula for the precession of the pericentre is derived. Then the gravita- 
tional radiation reaction is analyzed, where the famous Peters-Mathews formula 
is proved using the quaternion approach [7] . In this calculation we solved to can- 
cel the residual coordinate gauge freedom of the gravitational radiation reaction 
in the quaternionic equation of motion. 

The conclusion and the outlook is given in Section [5] followed by the Ap- 
pendix. 



2 Linearization and regularization with quater- 
nions 

2.1 Quaternion algebra basics 

The quaternion algebra has three imaginary units, generally called i, j and 6, 
where = = 4^ = — 1. Any of them anticommute 

ij = -jiJ« = -«j, H = -ie, (1) 

while the real unit 1 commutes with each of them. The four units together form 
the algebra's generators. Thus any clement of the algebra can be written afQ 

k = A + A^x + A2] + Ait = A + &. (2) 

Quaternion conjugation reverses the sign of the imaginary units 

A = A - a . (3) 

The magnitude of a quaternion is defined as 




Any quaternion can be written in the form A + aa, where a is the magnitude 
and a is the unit direction of a. Since a^ = — 1 the exponential on any quaternion 
can be expressed with Euler's formula 

gA+aa ^ gA ^^Qg a + a sin a) . (5) 

■^Note the unusual reversed ordering of the Ui components. 
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2.1.1 Representation of spatial vectors and rotations 

Spatial vectors are represented with pure quaternions, wliicli iias no real part 

yi = z\ + y] + xt, (6) 

where the z-axis is associated with the i unit. With this interpretation quater- 
nion multiplication can be expressed as 

AB = (A + a){B + b) = AS - a • b + Ab + Ba - a X b. (7) 

By substituting A = B = Q into this expression ([7]), the anticommutative cross 
product can be expressed as 

ab — ba , ^ 

a X b = . (8) 

Let us introduce a vector w. It is straightforward to show that a rotation around 
the vector w with magnitude |w| can be written as [4] 

i = RxM, (9) 

where R = e^. Note that MR = 1, hence x = RxR is the inverse rotation. Let 
us suppose that the rotation is parameter dependent R(s), and differentiate it 
with respect to s 

x' = RxR' + E'xR = xMR' + M.'M.x , (10) 



where definition ([5]) of x was applied. With the help of the identity (RR)' = 
R'R + RR' = and the cross product ([8]) this can be further written 

x' = iRR' - RR'i = Z X i, (11) 

where Z = 2RR'. Consequently Z is the instantaneous angular velocity of x with 
respect to s. With an inverse rotation the angular velocity Z can be expressed in 
a special coordinate system - in the Kepler frame ~ where x is instantaneously 
at rest 

Zk = RZR = 2R'R . (12) 
where the subscript indicates the Kepler frame. 

2.2 The Kepler problem with quaternions 

Equipped with the quaternion formulae we turn to regularize the Kepler problem 
with quaternions. The perturbed Kepler problem in the G = c = 1 system is 
given by the equation 

r+-3r = ef, (13) 
where r is the orbital separation vector of the orbiting bodies 

r = ri - r2, r = |r| , (14) 
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£ is the small parameter of the perturbation, and m = mi + m2 is the sum of 
the two masses. 

To regularize the Kepler problem the separation is defined by the following 
quaternionic equation 

r = UW, (15) 
where U = [/ + Usi + C/2j + J/ifi is a general quaternion. The conjugate of r is 

f = Uiu = UlU = -TOU = -r, (16) 

where the AB = BA, A = A and I = — I properties were used. Therefore r has 
no real part and can be written as r — zi + yj + xt. A direct calculation from 
tells us that 

x^U^ + U^- Ul - Ul , 
y = 2{UiU2 + U3U) , 

z = 2 {U1U3 ~ U2U) , (17) 

and the real part UUi — U3U2 + U2U3 — UiU = indeed vanishes. The obtained 
transformation p7p is just the KS transformation with the U ^ —U convention. 

m 

Transformation (|15p maps the four-dimensional quaternion space into the 
three-dimensional space of spatial vectors. Therefore the solution to a given 
three-dimensional r in terms of four-dimensional U is not unique. From (|15p it 
is clear that the transformation 

U ^ e'"U, (18) 

where a is an arbitrary real number, is a continuous symmetry of (|15p . It follows 
that we have a one-dimensional compact manifold - a fibre - of Us for a given 
r, and ([T8|) defines a fibration of the space of Us. The geometrical background 
of this transformation is elegantly described in Waldvogel (2005) [13]. This 
additional degree of freedom will be constrained in a careful manner. 

To complete the regularization the time has to be also transformed. The 
Sundman transformation is given by 

-^^2rJ-, 19 
ds \ m 

where s is the modified time and a is a real and at this point arbitrary function of 
s (it will be chosen such that it simplifies the solution) . From now the operator 
' indicates differentiation with respect to the modified time s. 

Inserting the definition of the orbital separation (1151) into the equation of 
motion (|13p while transforming the original time into the modified one using 
(fT9|) lead us to the following quaternionic differential equation [4] 

2U" - (2U'TF -ia)-+ 2W'- +m(-] - (v + W^] - + 4-eUrf = , 
r r V'"/ \ 2r J a m 

(20) 
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where 

r = mv'-ij'w = 2 {UiU' - uu[ + U2U3 - UsU^) , (21) 

which is a scalar in the sense that it is invariant under conjugation. This quantity 
is the four-dimensional scalar product of the tangent vector of the U (s) curve 
and the tangent of the fibre at that point multiplied by twclf|. Let us define a 
condition 

r = , (22) 

which means that the trajectory U (s) intersects the fibres under right angles. 
It is indeed a condition as transforming U with symmetry transformation (jl8p 
r transforms as 

r ^ r - 2aV, (23) 

hence condition can be satisfied by solving a first order differential equation 
for a{s). 

The geometrical constraint ([^ is equivalent with the so-called '^bilinear 
relation", which plays an essential role in the KS method[2]. It can be proved 
that (|18p is a dynamical symmetry, since the transformed U solves the equation 
of motion. Furthermore if condition p2|) is satisfied then T' also vanishes, 
which means that one can maintain this condition by finding the proper initial 
conditions U(0) and ILJ'(O) which satisfy the "bilinear relation" (|22p [i]. 

2.3 Solving the unperturbed case: Kepler orbits 

The unperturbed situation with condition (|22p reduces the perturbed equation 
of motion ([20)) to the following equation 

U" - (U'lJ' - 2a) - = . (24) 
r 

The coefficient of U is 

U'U' -2a _2a m\ _ 2ah 

r m \ 2 r J m ' 

where ft, is a constant of motion. Let us fix the parameter a 

which means that a is the semimajor axis of the elliptical motion. With this 
choice the equation of motion is the harmonic oscillator with constant frequency 

U[,' + Uo = 0, (27) 

where the subscript indicates the unperturbed case. 



This statement can be checked by differentiating the transformed U in JTSj with respect 
to a which produces the tangent of the fiber and then taking the four-dimensional scalar 
product with U'. 
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The general solution of this second order differential quaternion equation 
has six free parameters as it is constrained with (j22[) and has a redundant phase 
(fTB]). The trial solution of the unperturbed case can be parametrized as follows 

Uo=ai/2/3;i/'(q + /3g-i)M, (28) 

where 13± = 1 ± and q ^ e~ with uj = 2 {s — Sp) . In the next paragraph it 
is shown that the trial solution describes an elliptical Kepler orbit. 

Let us set R = 1 for the moment and substitute Uq into the definition of the 
separation 

ro = a/3+i«(z 4-/3^2-^ + 2/3) , (29) 

where z = and the identity qi = tq^^ was applied. This formula can be 
further expanded using the z = cos uj + i sin uj identity 

ro =a«(cosa; + 2/3/3^^) + aj;3_/3+^sinw, (30) 

which means that according to © ro describes the following parametric curve 

xq = a{cosuj + 2/3/3^""") = a(cosa; + e) , 

yo = a/3_/3^^ sincj = a\/l — sinoj , (31) 

with e = 2f3l3^^. These equations describe an ellipse in the {x,y) coordinate- 
plane, with semimajor axis a, and eccentricity e. From equations pip it follows 
that to = parametrizes the apocenter, thus w is equivalent with the eccentric 
anomaly, except that the latter is zero at the pericenter. This tells us that 
Sp is the time advance of apocenter passage measured in modified time. In 
the general case R is obviously the rotation between the orbital and reference 
frames, where the rotation according to ^ can be given with Euler angles 



(32) 



The formulae which provide the connection between the a, /3, Sp and angular 
parameters - the orbital elements - and the quaternion components are collected 
in the Appendix. 

Despite of the great advantages of the gauge condition (|22p - which is es- 
pecially fine for numerical calculations - for perturbative calculations another 
geometrical condition is proposed. 

2.4 The perturbed case 

The trial solution for the perturbed equation of motion (|20|) . is just the unper- 
turbed solution form (|28l) completed with general e proportional terms [4] 



U = aV2/3; V2 ( ^ ^ p^-i + + I ^ (33) 
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where both the ED and § quaternions are 0{e) quantities, and complex in the 
sense that they are in the subspace spanned by the units 1 and i, while the b 
quantity is real. 

Using the definition of the separation (jl5p the perturbed separation vector 
in the Kepler frame is 

vk = tkm + aP+^ {2t{z + - 2ilmE - 2ib} , (34) 

where the terms were neglected. The result tells us that parameter b describes 
a translation along the i unit which is a translation along the z-axis of the 
Kepler frame according to ([5]). By considering (|30p it parametrizes a translation 
perpendicular to the orbital plane. The same is true for the imaginary part of 
S, while the real part of S has no physical effect. Parameter B is complex and 
it is multiplied with t, thus the result is in the subspace spanned by the units 
j and 6. These units are associated ([6]) with the (x, y) coordinate plane of the 
Kepler frame, which is the orbital plane according to ([50]) . 

After introducing the trial solution in the perturbed case we fix the gauge. 
Vrbik's condition is that the real part of S must vanish[3] 

§* = -§, (35) 

where the operator * conjugates its complex quaternion argument. In this case 
the trial solution p3p has no t proportional part. 

Transformation (|18|) has a simple geometrical interpretation. It describes a 
double rotation, one in the (1,J) and another one in the (i,j) subspace. Hence a 
transformation (|18p on U ([55]) whose tangent is the coefficient of the t part of 
the trial solution (|33p divided by its real part cancels the coefficient of t. In the 
leading s order this rotation is 

(36) 



2(l + /3z-i)(l + /3z) ' 

which can be fulfilled without solving any differential equation for a in contrary 
to (USD. 

2.5 Example for the regularization 

The one-dimensional two-body problem is considered with the following special 
force / = ex^, therefore 

x+^=ex'^. (37) 

We have chosen this kind of special force since the equation of motion has a 
constant of motiorQ. 

* The first integral of the Eq. ^ is y'^ = Cy^e'^^y^ /4+m/2+emy'^e'^^y^ Ei{-2ey^), where 

OO 

C is the arbitrary constant of motion (C = 2Eq for unperturbed motion) and Ei(x) = — J 

— X 

t~^e~'^dt is the exponential integral function. 
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Eq. ([57)) can be regularized with the following transformations 

,2 



X 

dt 
ds 



y 



Then Eq. §7^1 is 



= 2eyiyr 



(38) 
(39) 

(40) 



In case of unperturbed motion (e = 0) the energy is the constant of motion 
{Eq = 2(?/')^/j/^ — m/y^) and Eq. (|40p is the equation of the harmonic oscillator 



0, 



(41) 



Cie''"' + C2e-'"^ where n = VT^blT^ is the orbital frequency. 



where clearly for bounded motion Eq < 0. The general solution for Eq. (|4ip is 



We assume that for perturbed motion (e 7^ 0) the form of the solution is 
Ue = Ho + sS. Then Eq. (|40)) in the leading order of (5 is 

S"-^-^S'-^S = 2y,{y',r, (42) 
2/0 2 

and using the yo ~ Acos{fls) solution of the unperturbed motion in (|42p we get 
S" + 2n tan{ns)S' + fl^6 = 2n^A^ cos(r2s) sin^ (fls), (43) 



where the sign of ft^S is positive, since Eq < 0. Numerical solutions for Eq. (|43p 
can be seen on Fig. [T] It can be seen that the numerical solution 6{s) of these 
two examples are well-behaving, bounded functions for various initial values. 




Figure 1: The homogeneous (left) and inhomogeneous (right) solutions for Eq. 
gal) for A = 1 = n, S'{0) = 0,(5(0) = 1 (dashed line) or J'(0) = 1,(5(0) = 
(line). 

So far we have represented the regularization of the one-dimensional per- 
turbed two-body problem with a heuristic special force. Let us consider the 
one-dimensional model using the generalized Sundman transformation 
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where x is not fixed yet. Eq. p7p can be regularized using transformations 
and (HI 



y 



iy' 



i'y' 



y X 2y° 



(45) 



li X = X (desingularized in perturbed orbit) we obtain Eq. (|42p . If x = xq 
{desingularized in unperturbed orbit) the result is 



y 



iy'T 



x'oy' 



262/ (y') 



Substituting the 



yo 



y xo 2y5 
-|- eSo perturbed solution, one obtains 

(^o oo = 2eyo (%) . 



(46) 



(47) 



where E{yo) ~ [2 (j/g)^ + 5m]/yg is not the constant of motion (note that the 
coefficient of the linear term is the constant of motion in case x = x (Eq. ([H])). 

In total two types of desingularization (i = x,xo) using the yo = Acos{ns) 
(we have fixed the frequency = 1) unperturbed solution can be given 



(5" + 2 tan(s)(5' + (5 = 2^^ cos s sin^ s, 



5m 



V 2^2 



sec s + tan s ] Sq 



2A^ cosssin^s. 



(48) 
(49) 



The numerical analyzis of these two equations with different initial values can 
be seen on Fig. 12.51 



5(s) 



<5(s) ' 



Figure 2: The numerical solutions for (|48p and with different initial values 

It can be seen that in this one-dimensional perturbed two-body problem the 
two types of desingularization methods lead to quite different solutions. There- 
fore the Sundman transformation is generally nontrivial in perturbed equations. 



3 Orbital elements with quaternions 

Before explaining the quaternion approach the classical equations are described 
in order to explain the relationship between the two different methods. The 
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equation system of the classical two-body problem is of total order six, hence it 
can be described with six first integrals, which are also called orbital elements. 
These elements are the semi-major axis a, the eccentricity e, the inclination 
9, longitude of the ascending node 4>, the argument of the pericenter ijj and 
the mean anomaly at the epoch Iq (or time of pericenter passage to) related 
to the dynamics. The Lagrange planetary equations in the standard perturbed 
two-body problem are|17) 

\ be smx + ^ — ^ — 



dt nVT- 
de VT^ 



2 



dt na 
dO r cos(x -I- ip) 



dt na^^/l- 
d(j} r sin(x -l- 



S* sin X + r (cos X + cos ^)] , 



dt na'^VT 

dip ^ d(h vT" 

= — cos ( ' 



T[l + ] sinx - S'cosx 



dt dt nae 

^ = -VT^^(^ + cose^)-s^, (50) 

dt \dt dt J na^ ^ ^ 

where x is the true anomaly, ^ is the eccentric anomaly and r is the parametriza- 
tion of the osculating orbit 

a(l — e'^) , , , , 

r = --^ ^ = al-ecosO, (51) 

1-hecosx ^ ^ ^ 

and I is the mean anomaly, which can be defined by the Kepler equation 

l -la = n{t- to) ^C-esin^, (52) 

and n ~ m^/^a~'^/^ is the mean motion. 

The S, T quantities arc the projections of the perturbing force to the orbital 
plane, while W is the projection to the normal vector of the orbital plane k 

5 = f-f, T = (kxf)-f, VF = k-f. (53) 

To derive quaternion differential equations for the orbital elements the trial 
solution p3p has to be substituted into the equation of motion ([20| . To sim- 
plify the calculation the quaternion equation of motion (|20[) is decoupled into 
two complex equations. The derivation of the complex equations are given in 
details |4j[T2 and the most important steps are briefly outlined in our Appendix. 
Here only the solution and the necessary definitions are presented. 



^in classical celestial mechanics the symbols arc f2 and uj respectively. Here we adopted 
the notations of J.Vrbik. 
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The following auxiliary quaternion quantities have to be defined 



^ ^ a Cx (rgfjf ) ^ a C'x (rgpfgo) ^ ^ / 2n ^ 

ml + /32 m 1 + Pz V/' 

W = -4e-rCx {ix) = -ie—roCx {ixo) + O (e^) , (54) 
m m 

where the operator Cx is a projector, which projects its quaternion argument 
to the complex subspaee spanned by the units 1 and i. The additional subscript 
indicates the unperturbed value of the symbol. The subscript K is omitted 
in r and rg as they are scalars and have the same value in every frame. To 
point out the relationship between the quaternion formulae and the classical 
equations (jSOp note that the real and imaginary part of the complex Q quantity 
is proportional to the previously introduced S and T (|53[) respectively, while W 
is proportional to W. 

The quaternion coefficients ([54| can be expanded into Laurent series 

n— -foo n— +00 

n— — ctt n— — oo 

together with D and S from ([33j) 

71— +00 n— +00 

B= ^ 2?„z", §= J2 (56) 

71— — oo n— 2 

The Laurent series are given in powers of z. From the definition of the orbital 
separation r (jisp follows that this is enough as the expansion of the separation 
contains every power of q. The coefficients and Si were left out 

from the expansion of D and § since they would only duplicate the q and q~^ 
terms of the solution ((33|) . while Sq was explicitly separated as b. 

Substituting expansions (|55|) and ((56)) into the complex equations (|92|) and 
([93]) the differential equations for the orbital elements can be extracted by 
matching the coefficients of z with the same power on both side of the equation. 



12 



The obtained differential equations are the following [T5] 



a' = 2aIm(Qo-/5Q-i), (57) 
13' = -^Im(Qi+3/?Qo + 3Q-i+/3Q_2), (58) 

Zi = -/3lilm(^^Wi+/3Wo) , (59) 

Z2 = -iRc(Wi), (60) 

Zs = ^Re{-/3+Qi+/3(l-3/32)Qo + (3-/3')Q-i 

+/3/3+Q_2}, (61) 
4 = f + ^Re{/3 (2 + /?^) Qi + (/?+ + 3/?4)Q„ 

-/3(l-2/?2)Q_i-/34g_^|^ (g2) 

5 = ilm {(/?2 Wo + 2/32^2) /3;i+/3>Vi}, (63) 
o 

and the two additional formulae for D and § is given in the Appendix. The 
Zi quantities are the components of Z/^, where the Kepler frame subscript was 
dropped to simplify the notation. Note that equation (|61|) is singular in /3, 
which shows that in the circular orbit limit the ordinary sense of the rotation 
no longer valid. 

The coefhcients Q„ or Wn of the Laurent series can be obtained with a 
contour integral, where Co is the unit circle 

Note that the Laurent expansion ((55)) of B and S has simplified the form 
of the differential equations (|57)) - ([63|) with respect to the Lagrange's planetary 
equations (|50p . The Laurent series of D and S absorbed the "short" term oscil- 
latory part of the equation. The remaining differential equations contain only 
the adiabatic, "long" term part, which might be easier to solve. 

We have to amend the equations above with the transformation of the an- 
gular velocity from the comoving Kepler frame to the inertial system, which are 
the following 

, Zi sin 'i/j + Z2 cos -ip 

<P = , 

smt' 

6 = Zi CDS')/; — ^2 sin?/!, 

ij}' = Z3-0'cos6i. (65) 

This transformation is familiar from classical mechanics, in deriving the Euler 
equations of the the rigid body. 
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4 GR perturbations 



In this section perturbations calculated from the general relativity arc examined 
using the described quaternion approach. The perturbations are examined with 
post-Newtonian approximation. The post-Newtonian approximation applies an 
expansion of corrections to the Newtonian gravitational theory with an expan- 
sion parameter e sa « m/r, which is supposed to be small, where v is the 
velocity. 

We use equations up to e^/^, (post)^/^-Newtonian order, which is the order 
where the dominant gravitational radiation damping forces occur. 

First of all, the (post) ^-Newtonian correction to the classical mechanics will 
be examined in the first section. This is followed by the (post)^/^-Newtonian 
analysis of gravitational radiation where we rederive the classical Peters-Mathews 
formula. 



4.1 Planar assumption 

The mentioned perturbations are planar perturbations, in the sense that the 
force lies within the orbital plane. In this case obviously S = 6 = and the trial 
solution (|33p contains only perturbations within the orbital plane 

Ua' = a^/"/?;'/' (q + I3q-^ + qB) , (66) 



It follows that in case of planar forces the §* = — § condition ([35|) is true. 
Remarkably the F = condition is also satisfied [3]. To show this we need the 
derivative of U expressed with Kepler frame quantities 

U' = V'k^ + UaM' = {v'k + Uk^^ (67) 



Therefore 

F = 2 Re (UW) = 2 Re jtUx^ [i^'k + ^k^^ k| , (68) 
and since the rotation can be dropped under the real part operator 

F = 2 Re (iJKiK + ^) ■ (69) 

In the planar case Vk is a complex number, therefore both Ua'Wj^ and 
tk are in the orbital plane spanned by the j and t units. In the planar case 
the orbital plane is preserved, therefore Zk must be perpendicular to this j, i 
subspace. It means that Z^- has only i part. It is easy to see from equation 
([69)) that the argument of the operator Re has no real part. Therefore in case 
of planar forces F vanishes. 

Consequently the equation of motion ([^0)) simplifies 

2U" - (2U'lF - 4a) - - U'- + 4-eUrf = . (70) 

r am 
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Let us introduce a r parameter by rescaling the modified time dr = 2a^/'^m~^^'^ds. 
With the help of the r parameter the equation of motion is just the perturbed 
harmonic oscillator ^ 

- hU + eUrf = . (71) 

In the planar case the equation of motion substantially simplified and identical 
with the equation of Waldvogel[T4]FI 

In the planar case the perturbations D can be expressed in a more conven- 
tional way. Let us introduce a = ao + 6a and /3 = /3o + <5/3 in the trial solution 
([28]) where 6a and 6(3 are first order quantities. In this case by matching the first 
order part of Eqs. (|28l) and ([M]) one obtaines the following important relations 

where A = l + /Sqz and B = z - /3o(l - f3o)'^A. 

4.2 The classical post-Newtonian effect 

In this section the leading contribution of general relativity to classical Newto- 
nian mechanics is examined in details. The force is given by numerous authors 

m 



apN = i n 



(1 -t- 3?7) - 2 (2 + ry) — - ^r/r^ 
r 2 



2 (2 - r]) fv } . (73) 



where the subscript PN denotes the post-Newtonian term, n — r/r, rj = 
{■mira2)lm? and v = |v| is the absolute value of the orbital velocity v =dr/dt. 
After transforming it to quaternion expression 

ap^ = - ^ J^(l + 3,7) (r^ri,) + ^2(2 + r,)m' 

^r.3^ ._^r^^ (74) 
8 a 2a 



This result ((74|) have to be substituted into the definition of Q ([54]) . where the 
separation tk , its magnitude r and their derivatives are treated as functions of 
z according to ((29)) . Therefore the result is a function of z 

+8zl3 (2 + p^Tj) + z2 {6 - 2/3^7] - 3) - 2?; + 0\32 + 677)}] . (75) 

Applying the contour integral (|64| the coefficients Q„ can be computed as 
follows. /3 < 1 therefore the only singularity of Q inside Co is at —[3. The other 
pole at —1/13 lies outside the unit circle. 



^Only the sign convention of h is different. 
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Q can be expanded around its pole at —(3 and keeping the coefficient of the 
(z + I3)~^ part, the resuh is 

Q_i = 2ma~^pz^l3+ {-^ - 8^^ - 3/3^ + 3/37? + 17/3^7; + f3^r]) . (76) 

In the same way with Q/z one finds 

Oo - -2ma~i/3Z^/3+ (-3 - 8/3^ - /?4 + ,^ + 17/32^ + 3/34^) ^ (77) 

Both of the coefficients Q_i and Qo are real. The differential equation for the 
semimajor axis is proportional to their imaginary part ([57)1 . therefore a' ~ 0. 
The remaining differential equations can be calculated in the same way. 

The resulting nontrivial differential equations in modified time for the orbital 
parameters are as follows [IB]. The equation for the argument of the pericenter 

6771 / /3_|_ ^ ^ 



*■ ^ ^ ^ I , (78) 



and for the modified time at apocenter 



= ^ (7/- 9 + /3^ (8?/- 15)) . (79) 



2a/3_ 



Using the transformation to the modified time (|19p in the leading order 

d 1 Pm d 
dt 2ay a ds ^ 

from (|78|) it follows that 

; 37773/2 
^= a^/2(l-e2) ' 

which is the known expression for the precession of the pericenter [20j. The 
remaining differential equations have zero on the right hand side of the equation 
and the corresponding orbital element is constant. 

4.3 Gravitational radiation reaction 

Gravitational radiation damping has been recognized as a process with very 
important observable consequences: the PSR 1913+16 system has given the 
first evidence that gravitational waves exist |21| , and other systems arc of high 
importance as well [22l [23] . The equation of motion is given by [24] 

8777772 

SiRR = ^ {-Ar,/^rh + S5/2V) , (81) 

07 

where the subscript RR indicates the radiation reaction term and 

A5/2 = 3(1 + + i(23 + 67 - 9/9)- - 5pr^ 
6 r 

B5/2 - (2 + 7)«' + (2 - 7)- - 3(1 + j)r\ (82) 



^ Our notation is slightly different from | 24| as the a and /3 parameters are occupied; 
instead we use 7 and p respectively. 
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The 7 and p parameters in represent the residue of gauge freedom that 
has not been fixed by the energy balance method and that has no physical 
meaning. It is known that these arbitrariness is equivalent with a coordinate 
transformation whose effect on the two-body separation vector is 



M5r = r + -f^[prr + (2p- 37)rv]. (83) 
15r^ 



We use transformation (jSSp to remove the gauge dependency from the quater- 
nion equation (PH)) . after substituting ([5T|) as the perturbing force. 

In order to apply transformation (jSSp on the quaternion equation of motion 
(|20[) we have to rewrite it in quaternion form using modified time (|19p. The 
definition of the modified time (|19p contains the separation r therefore any 
gauge dependent transformation of the separation, like leads to gauge 

dependent modified time s (7,p). Consequently the transformation of any real 
time derivative involves a new gauge dependent contribution 

d [m 1 d [rn 1 d pm I ^ d „ / c- 2\ /r, 
dt - H^rdS ^ H^rdS - H^-'^dS + ^ (''^ )■ 

It follows that transformation (|83|) in its original form docs not cancel these new 
gauge dependent contributions and needs to be reparametrized. The reparametrized 
transformation in quaternion form is the following 

^ Ua- + [K p r^Ug + {Lp-N^)rV'Kl (85) 

where K , L and N are unknown coefficients. 

To obtain them consider that Q ([54|) is Laurent series in z and gauge inde- 
pendence requires that every p and 7 proportional term in the coefficients must 
vanish. E.g. the coefficients of z^7, z^p and z'^0^p of Q ((54|) after simplification 
lead to a linear equation system which determines that TV = 3 and K = L = 1 
[18] . In the leading order according to (fT5)) this is equivalent with the following 
real time vectorial transformation 

S^lTTl^ 1 

r^r+^i-^[pfr+-(p-37)rv], (86) 

which is slightly different from (|83|) . 

The result from formula (|57 p -(|63 |) for the semimajor axis is now gauge 
independent [H] 

a' = T7T^ (6 + 97/32 _^ 219/34 + 97/36 _^ g^8^ ^ (37) 

15a'^'^ 



and also for the modified eccentricity 

15a5/2 



^ = 5/2 ('^6 + 273^2 ^ jg^g") . (88) 
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The remaining differential equations are trivial, with a zero on the right hand 
side, and the corresponding orbital elements remain constant. The gauge inde- 
pendent value of parameter D is given in [TH], while S is zero. 

Now we are in the position that the latter result for the semi major axis (|87p 
and also the expression for the eccentricity ((55)) can be easily verified. They 
must be equal with the two corresponding classical formula derived from the well 
known Peters-Mathews formula [7], which describes the effect of gravitational 
radiation. After substituting the expression e 

together with the transformation rule between the real and modified time (jit: 
one can derive the two equation below 

da 64 77to3 1 f , . 7^ 2 37 4 



dt 5 a3 (1 -£,2)7/2 \^ 24 96 

de 304 T^m^ e / 121 



^ + 777^ + T:^^ ] ' (89) 

(90) 



dt 15 a4 (1 _ g2)5/2 \^ 304 

which are indeed identical with the two formula derived from the Peters-Mathews 
equation [25]. 



5 Conclusion and outlook 

In this paper general relativity perturbations were examined using a new ap- 
proach where the regularization of the Kepler problem is given with quaternions. 
This approach is based on the usual Kustaanheimo-Sticfel method which is de- 
fined with matrices. 

With the new calculus the differential equations of the orbital parameters 
were derived in case when the perturbation is the leading (post)"'^-Newtonian 
order correction of general relativity. To test the new method the precession of 
the pericentre is rederived. 

Then the gravitational radiation reaction was analyzed, where the famous 
Peters-Mathews formula was reproved using the quaternion approach [7]. We 
have studied the gauge dependence of the equations of motion and we managed 
to remove the residual gauge freedom from the quaternionic equation of motion. 

The new quaternionic approach is easy to implement with program code. 
Quaternions can be represented with pairs of complex numbers, then the equa- 
tions can be calculated and solved with the help of complex analysis. This 
feature makes this method to a very efficient calculus for symbolic computa- 
tions. 

With the quaternion based regularization the spin-orbit and spin-spin inter- 
actions can be examined as well [5^. It is foreseen that these spin interaction 
related calculations provide the next step of our studies. 
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A The components of U expressed with orbital 
elements 

These formulae can be straightforwardly derived from the unperturbed solu- 
tion (|28p using the expression of the rotation ([5^ with the rotation angles 

U = cos ^ {cos {uj+ + I) + cos (uj+ - |) } , 

U,= a'/'p-'^'siJ-{cos[u^ + '^)+Pcos{u;^-^)} . (91) 
where w± = (0 ± i/") /2. 



B Decoupling the equation of motion 

For convenience the quaternion equation of motion (|20p can be decoupled into 
two complex equations. Premultiplying (^0]) with (—1 — U/f/ (2a) and also 
postmultiplying it by R while keeping only the 1, i part in 0{e) one obtains[l] 

- i (/?_ - /3z_) ^ + iz+/?' + 4i/3/3;i/3' 



+ (2/3_ + I3z_) Z3 - 4/3+s; + {1 + Pz)(b + 8z^ + 4^2^, 

^ \ az az'^ 

+ (1 + Pz-'^) B* + (1 - I3z) + 2z^^ + (1 - Pz-'^) + 2z 



dl 
dz 

(l + /?z-i)(l + /3z)2(Q, (92) 



where z± = z ± z ^ and Z„ are the components of the angular velocity vector 

m 

Similarly premultiplying equation (PH)) with (l + /3^) U^^/a and then keep- 
ing only the complex part in C(e) the second complex equation is the following 
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-8 



- 8 



z— + 8z— + 8z' 



+ /3z_(_ /3+ + dz dz 



1^ 



+ Z2Z^ 



/3/3+2+ + 2(l + /34) 



(1 + /3Z-1) (l + /3z)W(z). (93) 



C The solution for P and § 

Similarly by pairing the powers of z in the complex equations (j92|) and (1931) two 
additional equation for D and S can be gained 



1 



n— — oo 

{n + I) Q„ + (n + 1) 



n2 (n + 1) 



{n+l — 2 



1 «2' 



n (ri + 1) 



n2 (n+l)' 



'—Im 
4 



E 



(n— l)n — 1 n(n+l) 



z"(,94) 
(95) 



D The P and § quantity in case of the (post)^- 
newtonian effect 



The complicated quantity is given only up to second /3 order 

^(1 - 27?) + ^ (30 - 2z^ -9v + hzS) + 0(/?3) (96) 



while § is zero. 



E Gravitational radiation: P and § 

The fairly complicated D quantity is given in second j3 order 

V a 



16 ,^ /m\5/2 . 77/3^ /mN5/2 



+ (537 + 233z4) +0(/?3) 

45z^ V a / ^ ' 



(97) 



while § is zero. 
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